The ground state of a spin-crossover molecule calculated by diffusion Monte Carlo 
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Spin crossover molecules have recently emerged as a family of compounds potentially useful for 
implementing molecular spintronics devices. The calculations of the electronic properties of such 
molecules is a formidable theoretical challenge as one has to describe the spin ground state of a 
transition metal as the legand field changes. The problem is dominated by the interplay between 
strong electron correlation at the transition metal site and charge derealization over the ligands, and 
thus it fits into a class of problems where density functional theory may be inadequate. Furthermore, 
the crossover activity is extremely sensitive to environmental conditions, which are difficult to fully 
characterize. Here we discuss the phase transition of a prototypical spin crossover molecule as 
obtained with diffusion Monte Carlo simulations. We demonstrate that the ground state changes 
depending on whether the molecule is in the gas or in the solid phase. As our calculation provides 
a solid benchmark for the theory we then assess the performances of density functional theory. We 
find that the low spin state is always over-stabilized, not only by the (semi-)local functionals, but 
even by the most commonly used hybrids (such as B3LYP and PBEO) . We then propose that reliable 
results can be obtained by using hybrid functionals containing about 50% of exact-exchange. 



In nature there is a vast class of molecules whose mag- 
netic moment can be altered by an external stimulus. 
Typical examples of such molecules are the spin-crossover 
(SC) complexes [TJ [2], which, in their most abundant 
form, contain a Fe 2+ ion in octahedral coordination [3J 
and exhibit a transition from the low spin (LS) (singlet) 
ground state to a high spin (HS) (quintet) metastable 
state. Other examples are the cobalt dioxolene molecules 
[3HS] ■ These undergo the so-called valence tautomeric in- 
terconversion (VTI), namely an interconversion between 
two redox isomers, which differ in charge distribution and 
spin configuration. Both the SC transition and the VTI 
are usually observed for molecules in a single crystal and 
can be triggered by variations in temperature and pres- 
sure or by optical irradiation 7J . Furthermore it was also 
recently suggested that the VTI [5] and the spin ground 
state of a two-center polar molecule [9] can be controlled 
by a static electric field . 

SC complexes are promising materials candidates for 
molecular spintronics applications jTUl E] ■ Devices in- 
corporating such molecules are predicted to display dras- 
tic changes in the current-voltage curve across the phase 
transition |12[ 113] . and several transport experiments 
have recently achieved encouraging results. Alam et 
al. [2] were able to distinguish the spin state of a 
SC molecule placed on graphite by scanning tunnel mi- 
croscopy, while Prins et al. [15] demonstrated that the 
temperature dependent conductance of a device incorpo- 
rating a SC cluster correlates well with the phase tran- 
sition. In other cases, however, the data are not easy 
to interpret [16) and the experimental investigations are 
combined with density functional theory (DFT) simu- 
lations. In principle DFT should allow the computa- 
tion of quantities not easily accessible by experiments 
and should also provide parameters for effective trans- 



port models. However, unfortunately DFT results for 
SC molecules depend strongly (even qualitatively) on the 
choice of the exchange-correlation functional used [T71 - f22] 
and no standard has yet emerged. This essentially means 
that DFT is still not a predictive theory for this problem. 

Since most of the local and semilocal functionals un- 
derestimate the exchange energy, they tend to favor the 
LS state against the HS one [UJ [22] • This shortcome 
often leads to such large errors that even stable HS 
molecules are described as LS [THIUS]. In contrast, the 
most commonly used hybrid functionals are believed to 
over-stabilize the HS state [231 [23]. Besides, several au- 
thors [T5iri9, 22. 25 have criticized the common practice, 
that consists in assessing the performances of the various 
functionals by direct comparison with the experimental 
data. In fact, while experiments are usually performed 
for molecules in the condensed phase, DFT results refer 
to molecules in the gas phase. Since the properties of SC 
complexes depend strongly on environmental conditions 
(counter-ions, interaction between molecules, "strain ef- 
fects" etc...) [251 [27] their ground state may not be the 
same in these two phases. The question then becomes: 
can we produce a robust benchmark for DFT against the 
problem of predicting the Physics of SC compounds? In 
order to answer to this question ab-initio methods more 
accurate than DFT have to be considered. 

In the past wave-function based methods were used 
for this problem [16l 1X81421] . However, as the authors 
themselves pointed out, the results were plagued by sys- 
tematic errors ascribed to the limited basis set used for 
Fe 2+ and by the fact that the methods themselves ne- 
glect dynamic correlation (although this can be partly 
accounted for through a perturbative treatment). Here 
we propose an alternative route and perform diffusion 
Monte Carlo (DMC) [2H] calculations for a prototypical 
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Fe 2+ spin crossover molecule. As DMC represents one 
of the most accurate electronic structure method cur- 
rently available in order to compute ground state ener- 
gies, our calculations provided a solid benchmark for as- 
sessing the performances of DFT. In particular we con- 




FIG. 1: (Color on line) The cationic unit [FeL 2 ] 2+ (L=2,6- 
dypirazol-l-yl-4-hydroxymethylpyridine) used in the DMC 
calculations. Color code: C=yellow, 0=red (small sphere), 
Fe=red (large sphere), N=grey, H=blue. 

sider the molecule [FcL 2 ](BF4) 2 (L=2,6-dypirazol-l-yl-4- 
hydroxymethylpyridine) [25] (see Fig. [I]). We show that 
its ground state in the gas phase is HS but that a phase 
transition may exist in the solid state due to a number 
of crystal-related effects. We then show that the same 
result can be obtained by DFT hybrid functionals con- 
taining approximately 50% of exact exchange, thus con- 
firming early calculations for model molecules |25j . This 
establishes a recipe for the use of DFT for this class of 
materials, and it opens the opportunity to investigate 
with confidence the spin crossover transition of molecules 
in different environments (for instance on surfaces). 

DFT calculations are performed with the NWCHEM 
code [3U]- We use several functionals belonging to dif- 
ferent classes: 1) the Vosko-Wilk-Nussair local density 
approximation (LDA) [32]; 2) the generalized gradient 
approximation BP86, which combines the Becke88 ex- 
change functional [33] with the Predew86 correlation 
one [31]; 3) the hybrid functionals B3LYP [35], PBEO 
[351 137] and the Becke-HH [35], which include respec- 
tively 20%, 25% and 50% of exact exchange. We also 
consider a re-parametrization of the B3LYP functional, 
called B3LYP*, which includes only 15% of HF ex- 
change. This was introduced by Reiher and co-workes 
specifically in order to describe Fe 2+ complexes [33J [21] ■ 
The Ahlrichs triple-zeta polarized basis set [31] is used 
throughout. DMC calculations are performed by using 
the CASINO code [35J. The imaginary time evolution of 
the Schrddinger equation has been performed with the 
usual short time approximation and time-steps of 0.0125 
and 0.005 a.u. are used. Dirac-Fock pscudopotentials 
[40] |4T] with the "potential localization approximation" 
[412] have been used. The single-particle orbitals of the 
trial wave function are obtained through (LDA) DFT 
calculations performed with the plane-wave (PW) code 



QUANTUM ESPRESSO [13J. The same pseudopotentials 
used for the DMC calculations are employed. The PW 
cutoff is fixed at 300 Ry and the PW are re-expanded 
in terms of B-splines [H]. The B-spline grid spacing is 
a = 7r/G max , where G max is the length of the largest vec- 
tor employed in the PW calculations. Periodic boundary 
conditions are employed for the PW-DFT calculations 
and supercells as large as 40 A are considered. In con- 
trast, no periodic boundary conditions are imposed with 
DMC. Counter ions have been ignored and calculations 
are presented for the cation (Fig. [I]) to which we will 
generally refer as the molecule. 
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FIG. 2: Potential energy surface of the HS and LS state 
of a SC molecule. The collective coordinate r represents 
all of the nuclear coordinates of the molecule. The adi- 
abatic energy gap, A_E adla , and the vertical energy gaps, 
A£r S rt = A£ vcrt (r LS ) and A^s* = AS vert (r H s) are also 
indicated. 

The crucial quantity for understanding the spin 
crossover transition is the potential energy surface, 
schematically displayed in Fig. [2] This is typically plot- 
ted for the two different spin configurations as a function 
of a collective reaction coordinate r, which interpolates 
the molecule geometry along the LS to HS phase tran- 
sition. In our case DMC and DFT are used to compute 
the "adiabatic energy gap" 22 defined as 

A£ adia = E HS (r HS ) - E LS (r LS ) , (1) 

where r LS (r H s) and -Els(»"ls) [-^HS^hs)] represent re- 
spectively the geometry and the total energy of the LS- 
singlet (HS-quintet) state. When studying SC molecules 
at zero temperature Ai? adla is the central quantity, as 
it indicates whether the molecule ground state is LS 

( A£ ;adia > q) 0J . H g (A^^ia < Wg alsQ ca l cu l ate 

the "vertical energy gaps" [2"2~] 

AE™«(r LS ) = E HS (r hS ) - E LS (r LS ) , (2) 

A£ vcrt (r H s) = ^HS^HS) ~ ^LS(^HS) , (3) 

where -Ehs^ls) [-Els^hs)] is the energy of the quintet 
(singlet) state for the r^s (rns) geometry (see Fig. [2]). 
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Spin state 


d(A) 


BP86 


LS 


1.898(2), 1.958(4) 


BP86 


HS 


2.118, 2.127, 2.159, 2.166, 2.185(2) 


B3LYP* 


LS 


1.923(2), 1.996(4) 


B3LYP* 


HS 


2.182(2), 2.211(2), 2.207(2) 


B3LYP 


LS 


1.933(2), 2.011(4) 


B3LYP 


HS 


2.185(2), 2.218(2), 2.222(2) 


PBEO 


LS 


1.910(2), 1.984(4) 


PBEO 


HS 


2.169(2), 2.1965(4) 


Exp. 


LS 


1.909, 1.991, 1.912, 1.985, 1.980, 1.992 


Exp. 


HS 


2.105, 2.163, 2.103, 2.160, 2.170, 2.153 







A£™ rt (r L s) (eV) AE V 


er, (r HS ) 


eV) A£ adia (eV) 


DMC (At 


= 0.0125 a.u) 


0.65(3) 


-1.90(2) 


-0.36(4) 


DMC (At 


= 0.005 a.u) 


0.65(3) 


-1.91(3) 


-0.36(4) 



TABLE I: Experimental and calculated Fe-N bond-lengths 
for the [FeL2] 2+ cation. The number of bonds of a given 
length are indicated inside the bracket. The average difference 
between HS and LS Fe-N bond-lengths is about 0.2A, a typical 
value for SC molecules. 



DMC calculations were first carried out by using 
the molecular geometries optimized with DFT for the 
molecule in the gas phase. The Fe-ligand bond lengths 
computed with the various functionals arc listed in Tab.|T] 
As the DMC energy differences between the geometries 
calculated from BP86, B3LYP*, B3LYP and PBEO are 
of the same order of magnitude as the Monte Carlo sta- 
tistical error, we have not been able to firmly establish 
which functional produces the best structure. We have 
then decided to present results only for the structures 
relaxed with B3LYP, keeping in mind that the same are 
essentially valid also for BP86 and PBEO. 



Method 


A£J vert (r LS ) (eV) A£ v 


* rt (rns) 


(eV) A£ adia (eV) 


BP86 


2.87 


-0.180 


1.23 


B3LYP* 


1.97 


-1.013 


0.331 


B3LYP 


1.54 


-1.23 


0.012 


PBEO 


1.19 


-1.74 


-0.23 


Becke-HH 


0.51 


-2.50 


-1.33 


DMC (At = 0.005 a.u) 


0.28(4) 


-2.57(4) 


-1.19(4) 



TABLE II: Adiabatic and vertical energy gaps for the 
[FeL 2 ] 2+ cation calculated with DFT and DMC at the DFT- 
B3LYP relaxed geometry. The relative Monte Carlo statisti- 
cal error is indicated in bracket. 

The DMC adiabatic energy gap is reported in Tab. |TT] 
Our result indicates that the molecule in the gas phase 
is in its high-spin state, in contrast to the common belief 
and to the experimental result for the single crystal. Such 
ground state is indeed quite robust as DMC gives us an 
adiabatic energy gap of -1.20 eV. Since DMC provides a 
unequivocal assignment of the molecule ground state, it 
essentially establishes that no spin crossover transition is 
expected for [FeL 2 ] 2+ in the gas phase. Hence, in order 
to account for the experimentally observed SC transi- 
tion, one needs to understand how the embedding of the 
molecule in a crystal is able to reverse the relative order 
of the HS and the LS states at zero-temperature, i.e. to 
change the sign of Ai? adla . 



TABLE III: Adiabatic and vertical energy gaps for the 
[FeL2] 2+ cation calculated with DMC at the single crystal 
experimental geometry. We report DMC results for two val- 
ues of the imaginary time. The relative Monte Carlo statisti- 
cal error is indicated in bracket. The experimental molecular 
structure used for the calculation is taken from the X-ray data 
of reference [29] . 



The argument proceeds then as follows. Firstly, one 
has to repeat the calculations using the experimental 
geometries measured for the molecule in the crystal 
form |29) . These are less symmetric and present shorter 
metal-ligand bond-lengths than those optimized in the 
vacuum (see Tab. []]). The result is that the DMC- 
calculated A£? adla gets smaller, although it maintains 
the negative sign (compare Tab. Ill with Tab. III]). Sec- 



ondly, the electrostatic potential felt by the molecule in 
the crystal and due to both the counter-ions and the 
other molecules needs to be taken into account. This 
produces a relative shift of the HS and LS potential en- 
ergy surfaces. The magnitude of this effect, which tends 
to stabilize the LS state, has been recently estimated [STJ 
to be of the order of 0.5 eV. When the effect of the geom- 
etry and the electrostatic corrections are both included, 
DMC allows us to estimate a Ai5 adla for the condensed 
phase of about 0.2 eV. This is now positive, i.e. the 
ground state is LS, and very close to the typical values 
of the adiabatic energy gap inferred from experimental 
data [45] . 

We finally turn our attention to the assessment of the 
performaces of the various exchange-correlation function- 
als. Table [H] displays the vertical and adiabatic energy 
gaps calculated with DFT. We note that BP86 underes- 
timates the exchange so significantly that the molecule is 
predicted to be stable in the LS state ( A£ adia > 0) . Fur- 
thermore the absolute value of A.E vcrt (r LS ) [A.E vcrt (r HS )] 
is much larger than the corresponding one computed with 
DMC. This means that the standard local density ap- 
proximation predicts a very stable low spin ground state. 
B3LYP and PBEO improve only slightly the accuracy of 
the calculated gaps and the LS state still remains mas- 
sively over-stabilized. In contrast, as in the case of small 
Fe 2+ model complexes [25 , HH is found to be the func- 
tional that performs better, yielding a fair agreement 
with the DMC gaps. 

Importantly our analysis demonstrates that the assess- 
ment of the performances of a given DFT functional can 
be completely erroneous, if one insists on comparing the 
total energies calculated for the gas phase directly to 
experiments. If this was done with our DFT data, we 
would have concluded, as other authors did [531 HI]: that 
B3LYP* was the best functional. Our analysis instead 
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demonstrates that the correct assignment needs to be 
done against a reliable benchmark, with the result that 
the best suited functional must carry a fraction of exact 
exchange near to 50%. 

In conclusion, we have shown that for a typical SC 
molecule in the gas-phase the ground state is high spin 
and discussed how the this becomes singlet in the con- 
densed phase. We have also assessed the performaces of 
DFT against this problem, demonstrating that the HH 
hybrid functional, including 50% of exact exchange, is 
able to provide a quite accurate estimate of the energetic 
of the molecule. Our results then finally shed light on 
the long-standing issue of establishing the ground state 
of SC complexes. 

ACKNOWLEDGMENTS: We thank N. Baadji for use- 
ful discussions. This work is sponsored by EU (Hints 
project). Computational resources have been provided 
by TCHPC and ICHEC. 



[1] O. Kahn, Molecular magnetism (VCH, New York, 1993). 

[2] P. Giitlich, H.A. Goodwin, in Spin Crossover in Transi- 
tion Metal Compounds (eds. P. Giitlich, H.A. Goodwin) 
(Springer- Verlag, Berlin-Heidelberg, 2004). 

[3] SC molecules incorporating Fe 3+ , Co 2+ and Mn 3+ have 
also been reported. However these are much less common 
then those containing Fe 2+ ions. 

[4] O. Sato, J. Tao, Y.Z. Zhang, Angew. Chem. Int. Ed. 46, 
2152 (2007). 

[5] D.N. Hendrickson and C.G. Pierpont, Top. Curr. Chem. 

234, 63 (2004). 
[6] A. Bcni, C. Carbonera, A. Dei, J.-F. Letard, R. Righini, 

C. Sangregorio, L. Sorace, J. Braz. Chem. Soc. 17, 1522 

(2006). 

[7] A. Hauser, J. Chem. Phys. 94, 2741 (1991). 
[8] A. Droghetti, S. Sanvito, Phys. Rev. Lett. 107, 047201 
(2011) 

[9] N. Baadji, M. Piacenza, T. Tugsuz, F. Delia Sala, 
G. Maruccio and S. Sanvito, Nature Materials 8, 813 
(2009). 

[10] L. Bogani, W. Wernsdorfer, Nature Materials 7,179 
(2008). 

[11] S. Sanvito, Chem. Soc. Rev. 40, 3336 (2011). 

[12] N. Baadji and S. Sanvito, Phys. Rev. Lett., in press, also 

arXiv:1201.2028 (2012). 
[13] D. Aravena, E. Ruiz, J. Am. Chem. Soc. 134, 777 (2011). 
[14] M.S. Alam, M. Stocker, K. Gieb, P. Miiller, M. Haryono, 

K. Student and A. Grohmann, Angew. Chem. Int. Ed. 

49, 1159 (2010) 
[15] F. Prins, M. Monrabal-Capilla, E.A. Osorio, E. Coronado 

and H.S.J, van der Zant, Adv. Mater. 23,1545 (2011). 
[16] V. Meded, A. Bagrets, K. Fink, R. Chandrasekar, 

M. Ruben, F. Evers, A. Bernand-Mantel, J.S. Selden- 

thuis, A. Beukman and H.S.J, van der Zant, Phys. Rev. 

B, 83, 245415 (2011). 
[17] M. Swart, A. Groenhof, A.W. Ehlers and K. Lam- 

mertsma, J. Phys. Chem. A 108, 5479 (2004). 
[18] A. Fouqucau, S. Mer, M. Casida, L.M.L. Daku, 



A. Hauser, T. Mineva, F. Neese, J. Chem. Phys. 120, 
9473 (2004). 

[19] A. Fouqueau, M. Casida, L.M.L. Daku, A. Hauser and 

F. Neese, J. Chem. Phys. 122, 044110 (2005). 

[20] K. Pierloot and S. Vancoilie, J. Chem. Phys. 125, 124303 
(2006). 

[21] K. Pierloot and S. Vancoilie, J. Chem. Phys. 128, 034104 
(2006). 

[22] S. Zein, S.A. Borshch, P. Fleurat-Lessard, M.E. Casida 

and H. Chermette, J. Chem. Phys. 126, 014105 (2007). 
[23] M. Reiher, Inorg. Chem. 41, 6928 (2002). 
[24] M. Reiher, O. Salomon and B.A. Hess, Theor. Chem. 

Acc. 107, 48 (2001). 
[25] A. Droghetti, D. Alfe and S. Sanvito, in preparation. 
[26] M. Kepenekian, B. Le Guennic and V. Robert, J. Am. 

Chem. Soc. 131, 11498 (2009). 
[27] M. Kepenekian, B. Le Guennic and V. Robert, Phys. 

Rev. B 79, 094428 (2009). 
[28] W.M.C. Foulkes, L. Mitas, R.J. Needs and G. Rajagopal, 

Rev. Mod. Phys. 73, 33 (2001). 
[29] V.A. Money, J. Elhaik, M.A. Halcrow and 

J.A.K. Howard, Dalton Trans. 10, 1516 (2004). 
[30] M. Valiev, E.J. Bylaska, N. Govind, K. Kowalski, 

T.P. Straatsma, H.J.J, van Dam, D. Wang, J. Nieplocha, 

E. Apra, T.L. Windus and W.A. de Jong, Comput. Phys. 

Commun. 181, 1477 (2010). 
[31] A. Schafer, C. Huber and R. Ahlrichs, J. Chem. Phys. 

100, 5829 (1994). 
[32] S.J. Vosko, L. Wilk, M. Nusair and Can. J. Phys. 58, 

1200s (1980). 
[33] A.D.Becke, Phys. Rev. A 38, 3098 (1988). 
[34] J. P. Perdew, Phys. Rev B 33, 8822 (1986). 
[35] P.J. Stephens, F.J. Devlin, C.F. Chabalowski and 

M.J. Frisch, J. Phys. Chem. 98, 11623 (1994). 
[36] M. Ernzerhof and G.E. Scuseria, J. Chem. Phys. 110, 

5029 (1999). 

[37] C. Adamo and V. Barone, J. Chem. Phys. 110, 6158 
(1999). 

[38] A.D. Becke, J. Chem. Phys. 98, 5648 (1993). 

[39] R.J. Needs, M.D. Towler, N.D. Drummond and P. Lopez 

Rfos., J. Phys.: Condens. Matter 22, 023201 (2010). 
[40] J.R. Trail and R.J. Needs, J. Chem. Phys. 122, 174109 

(2005). 

[41] J.R. Trail and R.J. Needs, J. Chem. Phys. 122, 014112 
(2005). 

[42] L. Mitas, E.L. Shirley and D.M. Ceperley, J. Chem. Phys. 
95, 3467 (1991). 

[43] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, 
R. Car, C. Cavazzoni, D. Ceresoli, G.L. Chiarotti, 
M. Cococcioni, I. Dabo, A. Dal Corso, S. de Giron- 
coli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerst- 
mann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin- 
Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, 
A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, 

G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari 
and R.M. Wentzcovitch, J. Phys.: Condens. Matter 21, 
395502 (2009). 

[44] D. Alfe and M.J. Gillian, Phys. Rev. B 70, 161101(R) 
(2004). 

[45] G. Ganzenmiiller, N. Berka'inem A. Fouqueau, 
M.E. Casida and M. Reiher, J. Chem. Phys. 122, 
234321 (2005). 



